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Introduction 


This  paper  introduces  the  mixed  finite  element  method  as  a  viable  numerical 
procedure  for  the  boundary  contro’iability  of  the  linear  wave  equation.  Another 
numerical  implementation  using  Calerkin  finite  elements  has  been  investigated 
by  Glowinski,  Li,  and  Lions  in  [4].  However,  due  to  approximation  problems 
of  the  normal  derivative  on  the  boundary,  the  method  becomes  unstable  as  the 
mesh  is  refined.  To  correct  for  the  ill-posedness  of  the  approximate  problem,  a 
Tychonoff  regularization  method  was  implemented  in  [4].  The  aforementioned 
paper  also  presents  other  possible  remedies;  among  them  is  the  mixed  finite  ele¬ 
ment  method.  The  mixed  finite  element  approximation  is  a  plausible  procedure 
to  overcome  these  difficulties  since  the  derivative  at  certain  nodal  values  arises 
naturally  from  the  formulation. 

This  paper  is  numerical  in  nature;  related  theoretical  results  to  this  method 
will  be  presented  at  a  later  time.  The  first  section  gives  a  brief  description 
of  the  control  problem.  For  further  details,  we  refer  you  to  [4].  The  second 
section  of  the  paper  describes  the  mixed  finite  element  method  along  with  the 
approximating  spaces  used  in  the  procedure.  The  third  section  describes  how 
the  mixed  method  is  applied  to  the  controllability  of  the  wave  equation.  The 
last  section  presents  numerical  results  for  a  particular  test  problem  constructed 
in  such  a  fashion  so  that  the  exact  solution  is  known.  This  test  problem  was 
taken  from  [4]. 


2  Formulation  of  the  Control  Problem 


Let  be  a  bounded  domain  of  Rn  and  let  F  be  its  boundary.  Let  T  be  a  given 
positive  number,  where 


Q  =  x  (0,T),  E  =  T  x  (0,T). 


(1) 


Letp°eL2(«)  .p1  e  fi-’(fi). 

The  linear  wave  equation,  together  with  the  initial  conditions  p(x,  0)  =  p°(x) 
and  §f(x,0)  =  p1(x),  is 


(2) 


The  problem  is  the  following:  Is  it  possible  to  find  q  6  L2(E)  such  that 
adding  the  boundary  condition  p  =  q  on  E,  will  imply  p(x,T)  =  0  ,  §£(x,  T)  =  0, 
a.e.  ?  The  answer  is  yes  if  T  is  sufficiently  large.  The  proof  can  be  found  [5] 
and(6] . 

The  following  subsections  briefly  describe  a  method,  introduced  in  [5],  [6], 
and  [7],  for  constructing  a  boundary  function,  q  €  L2(E),  such  that  the  con- 
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ditions  previously  mentioned  hold.  This  again  has  already  been  presented  in 

[4]- 

2.1  Definition  of  the  Operator  A 
Define  E  by 

E  =  H^O)  x  L2(D);  (3) 

then  its  dual  E'  is  given  by 

E  =  H~\n)  x  L2(Cl).  (4) 

Now  define  A  €  L(E,  E  )  as  follows:  with 
e  =  (e°,  e1)  £  E,  solve  the  linear  wave  problem, 

ifitt  -  =  OinQ,  (5) 

V?  =  0on£,  (6) 

<P(*>  0)  =  e°(x),a.e.;  y>t(x,0)  =  e\x),a.e.  (7) 

Then  solve 

V’M  -  Arp  =  0  in  Q,  if  =  on  E,  ^(x,  T)  =  0,  a.e.;  rp,(x,  T)  =  0,  a.e.  (8) 
Finally,  define  A  by 

A  e  =  (0,(0),  —0(0)).  (9) 

The  fundamental  result  states  the  following: 

If  T  is  sufficiently  large,  then  A  is  an  isomorphism  from  E  onto  E‘ . 

The  proof  can  be  found  in  [6]  and  [7], 

2.2  Application  to  the  Boundary  Control  of  the  Wave 
Equation 

Let  £  £  E  be  defined  to  be 
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£=(p\-p0)-  (10) 

Now  consider  the  linear  problem, 

A  £  =  /•  (11) 

From  the  fundamental  result  it  follows  that  (11)  has  a  unique  solution  if  T 
is  sufficiently  large.  If  one  takes  the  solution  e  as  data  to  solve  (5)  —  (7)  and 

q  =  on  E  in  (8),  then  from  the  construction  of  A,  it  follows  that  p  =  ip  and 
p(T)  =  Pi(T)  =  0. 

It  can  be  shown  that  for  sufficiently  large  values  of  T,  A  is  strongly  elliptic 
from  E  onto  E  .  This  follows  from  [6]  and  [7J.  A  is  a  self-adjoint  operator, 
thereby  allowing  one  to  solve  the  problem  using  conjugate  gradient  methods. 
For  further  properties  of  A,  see  [4]. 


3  An  Explicit  Formulation  of  the  mixed  method 
for  the  linear  wave  equation 


Let  H(Cl;  div)  be  the  set  of  vector  functions  v  €  (L2(n))"  such  that  V  •  v  6 
L2(Cl).  Consider  the  linear  wave  equation: 

?£-Ap  =  0inQ  (12) 

with  p(x,0)  =  p°(x )  and  ff(z,0)  =  Pl(x).  Set  u  =  -Vp  .  Multiplying  by 
u  €  //(Q;  div)  and  integrating  by  parts  yields 


I  u  •  v  dx  —  I  pV  •  v  dx  —  —  I  p  v  •  rj  dy, 
Jn  Jn  Jr  " 


(13) 


where  g  is  defined  to  be  the  outer  normal  to  the  boundary  of  Cl.  Multiplying 
(12)  by  w  €  L2{Q)  and  integrating  gives 

J  ^jjwdx  +  J  V  •  u  wdx  =  0.  (14) 


The  system  is  then  approximated  using  finite  elements.  We  define  the  finite 
dimensional  subspaces  ,  V  and  W,  such  that  V  C  H(Cl\  div)  and  \V  C  L2(Cl). 
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For  convergence,  we  further  assume  the  property  that  ( V  ■  t;  |  v  e  V)  C  W. 

For  details  for  elliptic  partial  differential  equations;  see  [1],  [3],  and  [8].  For  the 
linear  wave  equation,  results  for  the  continuous  time  case  as  well  as  convergence 
results  and  stability  for  the  implicit  and  explicit  time  procedures  can  be  found 
in  Dupont,  Kinton,  and  Wheeler  [2].  In  this  paper  we  only  treat  the  following 
explicit  formulation. 

Spaces  satisfying  the  property  that  the  div  V  C  W  are  the  Raviart-Thomas 
spaces.  An  example  of  these  spaces  is: 


W  =  Mix(6r)  ®  Af}  i(5y)  (15) 

V  =  Mg(6t)®  Mi^6y)  x  Mlx{6x)®Ml{6y).  (16) 

Here  A/l_1(5r)®il/i.1(tfy)  is  the  tensor  product  of  piecewise  discontinuous  lin- 
ears.  Mq(8x)®  Af  lj($y)  is  the  tensor  product  of  piecewise  continuous  quadratics 
with  piecewice  discontinuous  linears.  Mi1(Sx)  <S>  Ml(8y)  is  the  tensor  product 
of  piecewise  discontinuous  linears  with  piecewise  continuous  quadratics. 

Let  At  >  0  and  tn  =  nAt.  Let  tt”  =  u(-,tn)  and  p"  =  p(-,  t"),  for  n  a 
positive  integer.  Define  (Un,  Pn)  €  V  x  Wby 


J  Un  ■  v  dx  —  J  Pn  V  •  v  dx  —  —  J  pn  v  ■  rj  dy,  V  v  E  V, 

r  pn  + 1  _  2 pn  i  pn- 1  r 

Jn - Af2 - wdx  +  y  y-fii)  =  0,Vu;ElV 


P1  -  P 
2  At 


P°  =  p(x,  0), 


-1 


(17) 

(18) 

(19) 

(20) 
(21) 


4  Discrete  formulation  of  the  Conjugate  Gra¬ 
dient  Method 


Recalling  from  Section  2,  e°  and  e1  are  initial  data  for  solving  the  forward  wave 
equation  (5)  —  (7)  so  that  the  normal  derivative  of  <p  on  E  is  the  boundary 
function,  q ,  such  that  p(T)  =  pi(T)  =  0.  In  this  discrete  formulation,  we  begin 
this  iterative  procedure  with  an  initial  guess  for  e°  and  e1.  The  subscripts 
denote  iteration  count. 


A 


Assume 


eg  6  W,  el0  £  W  (22) 

are  given; 

Now  solve  the  discrete  forward  wave  equation  : 


$n  ss  pn,  (23) 

Tn  fis  t",  where  rn  =  -Vp",  (24) 

f  T£  •  v  dx  -  /  $£  V  •  v  dx  =  0,  V  v  £  V,  (25) 
Jn  ~  ~  Jn 

r  *n+1  +  —  2*”  t 

/  -12 - - —wdx  +  /  V  •  TSu>dx  =  0,  Vit>6  W,  (26) 

At-*  Jn  ~ 

4>£+1  G  W,  (27) 

*o  !r=  0,  (28) 

n  =  0, 1 . N,  (29) 

where  the  forward  equation  is  initialized  by  $g  =  eo  >  ~  1  =  2 At  e£  -  Store 

<  ,<+1  ,  I?- 


Now  for  n  =  N,  N  -  1 . .  0  ,  compute  *J}  g  W  ,  T  |J  G  V  ,  £  W  by 

backward  time  integration. 

If  n  <  N  ,  compute  G  W  by  solving 


(30) 

(31) 

(32) 

If  n  =  N  ,  T  $  GVis  stored  from  forward  time  integration. 

Then  solve 


f 

Jn 


+  <*g+2  -  2<$g 


At 2 


’dx  +  f  V  •  Tg+1u>dz  =  0,  Vu;  £  W, 
Jn 

[  ■  v  dx  -  [  4>g  V  •  v  dx  =  0,  V  v  G  V, 

Jn  Jn 

$o  lr=  0. 
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/  ZS-  V  dx-  [  *2  V  •  v  dx  =  -  f  yon  V  •  2  dy,  V  v  £  V,  (33) 
Jn  Jn  Jr 

t  slf"-1  _i_  \Irn+1  _  Oxlr"  r 

— - — 175 - -wdx+  /  V  ■  ZSwdz  =  0,  Vwe  W,  (34) 

Jn  At  Jn 

K+1  ~K~1  =  *o  =  0,  (35) 

*o  lr=  Y0n,  (36) 

?o  «  £o.  where  £o  =  ~  Y^o>  (37) 

Xo  «  I  o.  where  I  o  =  ~  Y^o.  (38) 

*o"«-Io-2-  (39) 

Now  compute  g0  =  (0o>?o)  €  W  x  W  by  solving  the  discrete  Dirichlet 
problem, 


si  _  <*-i 


2o  Y  •  2  dx  =  0,  V  v  €  V, 

(40) 

wdx—  1  pxwdx,Vw£W, 

(41) 

Jn 

ffo  lr=  0; 

(42) 

and  then 


[  g^wdi  —  [  p°wdx—  f  dx,  V  w  6  W. 

Jn  Jn  Jn 


(43) 


If  g  o  =  0  or  small  then  set  e  h  =  e  o;  else  set  u>o  =  g  o- 
Then  for  k  >  0,  compute 

£  i+ii  £  ifc+i i 


(44) 


as  follows: 

'  Step  1:  Descent:  Tp  G  W  x  V. 
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[  T;  •  t ,  dx-  f  V  •  V  dx  =  0,  V  V  €  v, 
Jn~  Jn 


L 


*k  ±*k  ~2^t 

A  t2 


-r-0 


Store  <!>"  ,  +1  g  W  ,  T  ^  G  V. 

Now  for  n  =  A T,  TV  —  1 0  ,  compute  < 
backward  time  integration. 

If  n  <  TV  ,  compute  <^j!  6  W  by  solving 


> 

o' 

II 

€  V, 

(45) 

=  0,  Viu 

6 

(46) 

r=  0, 

(47) 

n  =  0,1, 

....TV, 

(48) 

(49) 

and  - 

-*Tl  = 

2  A  tw\. 

.  It  e  v 

\*T 

e  W  by 

L 


-r-n  ,  “r^+2 
*k  +  ~ 


At2 


wdx  + 


[  V  Tt+1 

•Tn  " 


u;  dx  =  0,  V  u;  6  W, 


(50) 


/  TJ-tnii-  /  $JV-  Wi  =  0,  VuGV  (51) 

'n  ~  Jn  ~  ~ 

$Tlr=0.  (52) 


If  n  =  TV  ,  T  ^  G  V  is  stored  from  forward  time  integration. 
Then  solve 
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J  ?k'vdx  —  J  'it  k  Y  ■  v  dx  =  —  J  Yl  v  ■  2  V  v  £  V, 


Jn 


-2»V 


wdx  + 


[  Y  Ik 

Jn 


w  dx  =  0,  V  xv  £  W, 


^  lr=  Y^l, 
Zk  «  ?J,  where?"  =  -  V^t, 

T*  w  I  J,  where  f  J  =  -  V^?, 

n^-iku- 


(53) 

(54) 

(55) 
1.56) 
(o7) 

(58) 

(59) 


Now  compute  gt  —  (StiSit)  €  x  W  by  solving  the  discrete  Dirichlet 
problem 


/  Q  t  ■  v  dx  -  f  jj  V  •  v  di  =  0,  V  »  £  V, 
Jo  Jn 

/  V  ■  Qkwdx  =  f  ^k-.^k  wdx,Vw£  W, 
Jn  Jn  *At 

^  lr=  0; 


and  then 


/  </*u>da:  =  - 
Jn  Jn 


—  j  &  *  U/  dx ,  V  Ztf  £  l  K 
n 


Then  compute  by 


(60) 

(61) 

(62) 


ffiO\ 
\  ■>°/ 


(64) 


/>*  = 


fn  9*  '  ^  kdx  +  fnglkgldx 


In  §*  •Ii(/z  +  fnJkwkdx 
Once  /j*  is  known,  compute 


(65) 
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(66) 


e  t+i 

$t+i 
«i+i 

£  k  +  \ 

If  a  *+i  =  0  ,  or  is  small  ,  then  set  e  h  =  e  t+i  ,  =  ^t+i,  =  'J't+i  i 

else  compute 


=  e  t  ~  pk  wk, 

=  $k-Pk$k,  (67) 

=  tfi-p*#*,  (68) 

=  £k-pk£k-  (69) 


fn  ©k+i  •  ©k+i<tr  +  /nSjt+1si+ici* 

7k  =  - 

fndk-  0 kdx  +  Jn g\g\dx 


(70) 


Set 


U>  fc+i  =  a  t+i  +7t  \Vk  (71) 

and  k  =  k  +  I  and  go  to  Step  1. 

Remarks: 

•  As  pointed  out  in  [4],  substantial  computer  memory  cost  is  reduced  by 
solving  the  wave  equation  backward  in  time. 

•  This  formulation  is  only  valid  for  problems  with  smooth  data.  A  variant 
of  this  conjugate  gradient  method  is  required  to  handle  nonsmooth  data. 
Procedure  can  be  generalized  to  treat  this  case. 


5  Numerical  Results 

This  method  was  used  for  a  test  problem  used  in  [4].  In  [4],  an  exact  solution 
is  constructed  for  the  problem  A  e  —  £  on  the  unit  square.  For  details  of  this 
calculation,  we  refer  you  to  [4].  Only  the  results  will  be  presented  here. 

If 

e°(j)  =  sin  sin  WX2,  (72) 

e‘(x)  =  7r\/2sin  ?rxi  sin  n-x;,  (73) 
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then 


ip(x,  t)  =  \/2  cos  irv/2 (t  - 


1 

4y/2 


)  sin  irxi  sin  tt:e2. 


Defining  r,  ,*  =  1,5, 3, 4  by 


r,  =(i|ie  r,xi  =  o), 

r2  =  (x  |  x  e  r,zi  =  i), 
r3  =  (x  |  x  G  r,z2  =  0); 
r4  =  (x  |  x  e  r,ar2  =  i), 


(74) 


we  have 


|r,ur3=  — t\/2  cos  ?r\/2(i  -  ^=)  sin  7ra:2, 
^  |r*ur«  =  -7r\/2cosir\/2(<  -  -^=)sinjrx!. 


(75) 

(76) 


Using  final  time  T  -  ^(n  +  |)  (n  is  a  nonnegative  integer),  ip  =  ip0  +  xpi , 
where  0O  and  0i  are  the  following: 

0o  =  -7T\/2cos  jr\/2(t  -  ^—^(sinirxj  cos27rx2  +  cos2irxlSin7rx2),  (77) 

0!  =  47r(T  —  t)  sin  7t\/2(<  —  •—=)  +  (78) 

28 

(_  1)"  3\/I  S*n  -  ^)s*n7rxi  sin  7rx2 


+4sinirxi  ^ 


m 


m  >3 

m  odd 


m2  —  1 


2(— l)n+1 


r— — ==-sinir\/ 1  -f  m2(t  -  T) 
v  1  +  m2 


3V2 

H - z - 7  cos 


2  i  vv‘'  —  - —  ) 

m2-4  4>/2 


sin  mnx2 


+4sin7rxo  ^ 


m 


m  >  3 

m  o«id 


m2  —  1 


-1V+1 


2( —  1 ) 


l  vT+  m2 


Sin  7T  \/l  +  m2(<  -  T) 


+ 


3\/2 
m2  —  4 


>*V2(t  -  -L) 

4v/2 


sin  rmrx\. 


Since  p  =  0,  we  compute  p°  and  p1  from  0O  and  0j  so  that 
P°(x)  =  0o(-c,O)  +  0i(x,O), 


p'W  =  **<x,0)+<2L(x.„). 


501 


(79) 

(SO; 
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Since  p°  and  p1  involve  infinite  trigonometric  series,  Fast  Fourier  Transforms 
are  used  for  these  calculations  (  m  is  taken  to  be  255). 

In  the  conjugate  gradient  algorithm,  e°  and  e1  are  initialized  to  be  ze>-o  and 
final  T  is  (n  =  3).  The  following  pages  represent  calculations  for  h  =  1/16 
,  1/32,  and  1/64.  The  first  six  plots  represent  graphs  cf  the  calculated  e°  and 
e\  along  with  the  known  e°  and  e1.  The  last  three  plots  represent  variations  of 
I!  ?(0  Ili’Cr)  and  II  £c(0  ||z,J(r)  with  t.  All  approximate  solutions  are  represented 
by  dotted  lines  and  known  solutions  are  represented  by  solid  lines. 


The  first  table  shows  that  the  method  is  much  better  behaved  as  the  mesh 
is  refined  However,  the  second  table  shows  that  the  iteration  count  goes  up  as 
the  mesh  refined;  roughly  speaking  like  h~  5.  This  is  substantially  better  than 
the  Galerkin  finite  element  procedure  without  the  regularization  discussed  in 

[4]- 


h  = 


i/8 _ _ 

3.03  x  10~2 
5.69  x  10-" 


1/16 

1.00  X  10-‘j 


1/32 


1/64 


~  e°c  |U^(n) 
1  ‘ka(n) 


3.11  x  10~3 


1.25  x  10 


eL  —  e 
e°  —  e 


1.38  x  lO"1 
2.85  x  10- 2 


1.79  x  10'^ 
4~95  x  10-5 


9.76  x  10~3 


4.22  x  10~3 


<lc  IIl^S) 


1.02  x  'l0“f 


1.70  x  10-2 
3.31  x  10~3 


7.39  x__10~3 
1.37  x  10~3 


h 

no.  of  iterations 

1 

X 

19 

i 

a 

30 

16 

48 

1 

-32 

72 

1 

64 

119 

7.298 


7.401 


7.394 


6  Conclusion 

From  the  numerical  results,  the  ill-posedness  of  the  approximate  problem  is 
allieviated  considerably  when  using  the  mixed  finite  element  procedure.  Even 
though  the  iteration  count  goes  up  as  the  mesh  is  refined,  there  is  no  oscillatory 
behavior  present  as  in  [4]  with  no  regularization. 
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1/16 


=  1/32 


1/32 


g(t) |  L2  on  boundary  where  h  =  1/16 


♦  * 
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iteratively  adjusts  the  matching  conditions  at  the  interfaces  of  the  subdomains.  Numerical  results  are 
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0.  Introduction.  In  (lj  Glowinski  and  Wheeler  defined  domain  decomposition  algorithms  for  solving 
mixed  finite  element  approximations  of  elliptic  problems  with  non-constant  coefficients.  A  key  result 
in  [1]  was  the  formulation  of  the  matching  conditions  at  the  interfaces  of  the  subdomains  as  variational 
problems  defined  over  convenient  trace  space.  These  new  problems  were  solved  by  conjugate  gradient 

_  .  C 

algorithms  using  simple  preconditioners  resulting  in  a  0(h  ’  )  number  of  iterations  to  achieve 

convergence.  In  this  paper  we  shall  discuss  a  procedure  for  accelerating  the  convergence  of  the  abvA'e 
algorithms  which  is  essentially  based  on  a  multi-level  technique  acting  on  the  trace  space  associated  to 
the  interfaces. 

In  Section  1,  we  shall  give  some  examples  of  elliptic  problems  originating  from  flow  in  porous 
media.  Compared  to  more  traditional  solution  methods  the  algorithm  described  in  this  paper  have  been 
quite  successful  as  we  shall  demonstrate  in  Section  4.  In  Section  2  which  follows  closely  [1]  we  shall 
recall  the  mixed  variational  formulation  of  elliptic  problems,  the  mixed  finite  element  approximations 
and  the  associated  domain  decomposition  methods.  In  Section  3  we  shall  discuss  a  multilevel  method 
to  speed  up  convergence  of  the  domain  decomposition  algorithms  discussed  in.  Section  2.  Results  of 
numerical  experiments  will  be  discussed  in  Section  4.  Finally  some  mesh  refinement  methods  w’ell 
suited  for  domain  decomposition  and  mixed  finite  element  methods  will  be  discussed  in  Section  5. 

1 

1.  Motivation  for  Robust  Elliptic  Solvers. 

In  our  first  example  we  consider  the  pressure  equation  which  arises  from  miscible  displacements 
in  porous  media.  The  equation  has  the  form 

(1.1)  u  =  — A  grad  p  in  Q, 

(1.2)  V-u  =q  in  fl, 

(1.3)  u  •  v—  0  on  dQ, 

where 

A  =  k(x,  y)/p(c). 

0 


In  this  problem  fl  is  the  flow  region,  u  is  the  Darcy  velocity,  p  is  the  pressure,  q  is  a  source  or  sinks 
term,  k  is  the  permeability  of  the  porous  media,  p  is  the  viscosity  of  the  concentration  c  of  the  fluid 
which  is  flowing  through  the  porous  media.  In  this  example  we  use  a  permeability  field  and  a  form  of 
the  viscosity  which  has  been  previously  obtained  from  laboratory  experiments.  In  Figure  1.1,  a 
visualization  of  A  is  shown.  In  this  case  we  have 

min  A  =  .810xl0  —  “  and  max  A  =  . 282x10“^, 

implying  that  (1.1)-(1.3)  is  badly  conditioned.  However,  as  it  will  be  seen  with  more  detail  in  Section 
4,  we  have  been  able  to  solve  this  problem,  using  domain  decomposition,  in  less  than  10  iterations. 


Variation  of  coefficient  A 
Figure  1.1 


3 


2.  Mixed  Formulation  of  Elliptic  Problems  •  Associated  Finite  Element  Approximation  and  Domain 


Decomposition. 

2.1  The  Model  Problem. 


We  consider  on  ficRn  the  following  Neumann  problem 
f— V-AVp  =  f  in  Q, 


(2.1) 


[AVp-t/  =  g  on  dCl(  =  T), 

where  v  is  the  outward  normal  vector.  We  assume  the  compatibility  condition 


(2.2)  J  fdx  +  j  gdF  =  0. 
ci  r 

Our  formalism  is  motivated  from  How  in  porous  media  where  (2.1)  is  the  pressure  equation,  but  the 
method  to  be  described  applies  to  other  branches  of  science  and  engineering.  Also  we  have  been 
considering  the  pure  Neumann  problem  since  it  is  the  one  occurring  most  frequently  in  applications.  In 
fact,  it  is  also  the  most  difficult  case. 

2.2  A  Mixed  Variational  Formulation  of  Problem  (2.1) 

Define  u  by 


(2.3)  u=  —  AVp. 


We  then  have 


(2.4)  V-u-f=0, 


and 
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(2.8)  V0  =  {v  |  v«H(n,  div),  v<i/  =  0  on  T). 

Here 

(2.9)  H(D;  div  )  =  {  v  e(L2(n))n  and  div  v£L2(Q)}. 

_1 

Suppose  fcL2(n),  ge  H  2(r)  and  A  is  symmetric  such  that  Ae(L°°(n))n:<n  and 

A(x)£-$>a|$!2,  V£eRn,  a.  e.  on  o 

with  a  a  positive  constant. 

If  (2.2)  holds  then  (2.1)  has  a  unique  solution  on  H^(fi)/R  implying  the  uniqueness  of  u.  An 
alternative  formulation  of  (2.1)  is  given  by 

Find  pfL2(fi),  ucH(Q;  div),  such  that 
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Examples  of  particular  finite  element  spaces  for  which  (2.11)  is  well  posed  and  for  which 

*im  uh~ *u  an^  ^;rn  P'n-"1 ‘P  can  'Be  f°und  in  [2].  Additional  convergence  results  including  error  estimates 
can  be  found  in  [3,  4]. 

2-4  Domain  Decomposition  Method  for  Problem  (2.1),  (2.11). 

IVe  follow  here  the  notation  and  methodology  developed  in  [1],  Considering  first  the 
continuous  problem  whose  tormula  is  much  simpler  we  suppose  that  Cl  has  been  decomposed  in  two 
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subdomains  and  fl0.  Figures  2.1a  and  2.1b  show  such  domain  decompositions  and  corresponding 
notation. 


If  we  denote  by  {p^,  u-}  the  restriction  of  {p,  u}  to  Q.  there  is  equivalence  between  (2.10)  and 

J  (V-m— f)  qjdx  =  0,  V  qjcL2(a), 
a 

(2.12)  .  i 

1  /(A-1ui.vi-PiV.v.)dx  =  0)  Vv.eVio,  i  =  l,  2, 

•a 

1 

(2.13)  Uj-iaH hg  =  0  on  FnSQ-,  i  =  1 ,  2, 

(2.14)  J2  ui  =0  on  T> 
i  =  l 

(2.15)  £  J  dx  =  0,  VvcV0, 

1  -  1  n 


l 


Vio  =  {vjlvifII(^j>  d»v).  Vj-J/j=0  on  dfij. 

Since  V0  =  V^0©V20©V^0  (where  Vy0  is  a  complementary  subspace  of  Vjo®  \;20  in  VG)  it  follows 
from  (2.12)  and  (2.15)  that  (2.15)  can  be  replaced  by 

(2.1C)  E  /  (A-1  uj-v-pj  V-v)  dx  =  0,  VvcV70. 

',  =  1  fii 

In  addition  to  (2.12)-(2.1G),  {p-,  m}  must  satisfy  the  compatibility  condition 


(2.17)  J  fdx  +  j  g  dT  +  J  Uj-j/.  d7-  =  0. 


a  5a  nr 


From  (2.12)-(2.16),  the  local  solutions  satisfy  at  the  interface  7  the  matching  conditions  (2.14) 

and  (2.16).  From  this  observation  we  can  generate  two  classes  (at  least)  of  iterative  methods  for 

solving  problem  (2.11)  by  domain  decomposition.  In  both  approaches  we  assume  that  one  of  the 

matching  conditions  is  satisfied  by  an  appropriate  choice  of  boundary  conditions  on  7  and  we  try 

iteratively  to  satisfy  the  other  matching  condition.  In  this  paper  we  shall  concentrate  on  the  case 

where  the  balance  given  by  (2.14)  is  satisfied;  we  try  therefore  to  verify  (2.16). 

1 

This  leads  to  the  introduction  of  a  variational  problem  involving  functional  spaces  defined  on 
7.  Precisely  such  a  functional  space  is  V~0  defined  by 


(2.18)  V°0  =  {fi\v  (\'7o,  J  fi-v  d7  =  0). 


We  define  next  a  bilinear  form  a(-,-)  over  V°0  x  V°0  as  follows: 

Consider  ^rV°0;  we  associate  to  fi,  u ■  (/j)  and  Pj(p)  by  solving 
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(2.19) 


J  V.ujM  Vjdx=0,  VvifL2(fii), 

°i 

(2.20)  J  (A.-1  Uj(/i)  -Vj— pj(p)  V-Vj)dx=0,  Vvi£Vio, 

fii 

(2.21)  Uj(p)  -Vj=0  on  rn3flj,  Uj(/i)  •  lAS/i-i/.  on  7. 

Since  /  u-(/i)-i/-  dF- =0,  the  above  problem  is  well  posed  in  11(0-,  div)  x  L2(fb)/R.  To  insure 

dfij 

uniqueness  of  Pj(p)  we  enforce  the  conditions 

(2.22)  y  p1(^)dx=0,  £  y  (A'1  Uj^-n-Pjtp)  V-n)  dx=0 
«1  i  =  1 

where  IIe(V^0  —  V°Q).  Finally  we  define  a(v)  by 

(2.23)  a(/i,  p')  =  £  J  (frr'-Piil*)  V-/)dx,  V  ,/£V°0. 

1-1  ni 

It  has  been  shown  in  [1]  that  the  bilinear  form  a(-  ,  •)  is  symmetric  and  positive  semi-definite  over 
V°0xV°0.  Moreover,  it  is  elliptic  for  the  norm  induced  by  H(fl;  div)  over  the  quotient  space  V°0/R, 

A  A  .  / 

where  R  is  the  equivalence  relation  defined  by  fj.  R  n  [fx—^x  )-i/  =  0  on  7. 

From  the  above  result  we  can  interpret  (2. 1 2)-(2. 17)  as  a  linear  variational  problem  in  V°0. 
To  formulate  this  latter  problem  consider  A0  £H(Q;  div)  such  that 

(2.24)  Ao-i/  +  g  =  0  on  F, 

(2.25)  |fd*-F  /  gdr+  j  Ao.i/idT  =  0,  Vi  =  l,  2; 

rnsa  7 
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solve  then  for  i  =  l,  2, 


(2.2C)  /(V.uoi-Oqid>:=0,  VqifL2(I)i), 

(2.27)  J  (A"1  uoi-vi-poiV.vi)dx  =  0,  V  VjtV^, 


(2.28)  uoi-,/i +e=°  °n  7(1^, 


(2.29)  uoi  - ^  =  A0  - on  y. 

The  constants  associated  to  the  pQj  are  adjusted  as  follows: 


(2.30)  J  poldx=0, 

% 

(2.31)  £  [  (A-1uoi.n-poi  V-n)  dx=( 
i  =  l  A 


Let  us  now  denote  by  u0  the  element  of  H(fJ;  div)  such  that  u0|q  =:u0j-  If  we  define  u  by 


(2.32) 


we  cleariy  have  ucV0.  Denoting  Ac\%0  as  the  component  of  u  in  the  decompositic 
V0=Vj0©V20©V-r0  we  have  from  (2.17),  (2.25),  (2.28),  (2.29)  that 


(2.33) 


A-ia  d7  =  0,  i.t.  AcY°0; 
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define  similarly  pj  by  p-  =Pj  —  P0j  • 

We  have  then 

(2.34)  J  V-uiq;dx  =  0,  V  q;fL2(a), 
a 

(2.35)  J  (A-1  uj-vj-pj  V.Vj)dx  =  0,  VvjcV^, 
fii 

(2.36)  iij  •  ia=0  on  dfbnr,  u-  -ia  =  A-tA  on  7, 

(A  —  1Uj’IT— f>.  V  n)dx  =  0. 

It  follows  from  (2.16)  that 

(2-3S)  t  I  (A_1u;  •  ,.-p,  V.,,)  dp  =0,  V  p  (V°0. 

i=14 

From  the  definition  of  u^,  p-  and  from  (2.38)  we  obtain 

(2.39)  £  /( A-1nj-/i-j5j  V-p)  dx=-£  /  (A_1uoi-/j-poi  V-/i)  dx,  V//V°0. 

i=1  n.  i=1  a 

It  follows  from  (2.23)  and  (2.33)  that  A  is  the  unique  solution  of  the  linear  variational  equation 

,  Find  AeV°0  such  that 

(2.40)  ■ 

a(A,  p)  =-£  /  (A-1uoi-Ai-poi  V-M)dx,  V  /rtV%. 
i  =  1  O. 


(2.37) 


i.x =0,  £  t 

i  =  1  A 
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In  [I],  we  showed  that  the  variational  problem  (2. 40)  can  be  approximated  by  a  finite 
dimensional  problem  of  the  same  nature,  obtained  by  combining  the  mixed  approximation  of  Section 
2.3  with  the  domain  decomposition  principle  of  Section  2.4.  In  addition,  a  conjugate  gradient  method 
for  solving  this  finite  dimensional  problem  approximating  (2.40)  was  discussed  in  detail  in  the  above 
reference. 

In  the  following  Section  3,  we  shall  describe  multilevel  techniques  for  solving  the  finite 
dimensional  problem  approximating  (2.40);  it  can  be  seen  as  a  multigrid  method  operating  on  the 
interface  7. 

3.  Multilevel  Solution  of  Problem  (2.40). 

3.1.  Domain  Decomposition  of  the  Discrete  Problem. 

Following  Section  2.3,  it  is  easily  shown  that  the  discrete  mixed  problem  (2.11)  is  equivalent  to  finding 


K,i’ 

ph,i } >  i  =  l>  2,  satisfying 

(3.1) 

/<v' %,i-0 <^=°.  vvwh., 
a 

1 

(3.2) 

/(A_lu h,i*  vi  -  Ph,iV'vi)  d*  =  0, 

n. 

1 

VVVo  h,i 

(3.3) 

/  K.i-^+s)  v-i/ dr=o,  Vvitvoh., 

aanr 

(3.4) 

0 

E  uh  i-I/i=0  on  7' 

i  =  l  1 

(3.5)‘ 

.E  /  (A_luh,i'v-Ph,i  v-v) dx=0- 

VvcVoh- 

n 


where  VQjJ  j  (resp.  j)  is  equal  to  V  ,  U  (resp.  W.  Jq  ).  As  in  the  continuous  case  we  associate  to 
’  :  i  i 

7  a  complementary  subspace  ^  of  0  in  V^;  that  is 
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Voh  ~  Voh,l  ®  Voh,2®  Voh,7‘ 


It  follows  from  (3.1)  and  (3.2)  that  (3.5)  can  be  replaced  by 


(3-6)  t  J  (A  ^h/v-pb-  V-v)  dx=0,  V  vcVoh)7. 

1-1  n. 

1 

In  addition  to  (3.5)  and  (3.G)  {u^  j,  p^  j}  has  to  satisfy  the  compatibility  conditions 

(3.7)  J f  dx  +  J  gdT  +  J u)1j-";d7  =  0,  i  =  l,  2. 
a  5a  nr  r 


Finally  we  decompose  ^ 


as  the  direct  sum, 


(3.8) 


V  =  V° 

oh, 7  oh, 7 


®  V 


n 

oh, 7 


where 


(?••)  VSh,T  =  (“V„h,T  I 


(3.10  )  ^  o^h  =  f  ^  and  n<r  VQb  ^.wiih  J  U-u  d 7 9^ 0} . 

7 

3.2.  Discretization  of  the  Boundary  Problem  (2.40). 


Following  the  development  in  Section  2.4,  we  approximate  (2,40)  by  the  following  variational  problem 
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in  V°L  „  x  V°v 


oh, 7  oh, 7" 


(3.11) 


Find  A^  <  V®h  such  ihai 


ah  =  J  (A  \Wl  ■  /i-Pohii  ^-/Odx,  V/i«V° 

1  =  1a 


where  A^,  u^  j  and  p^  j  are  obtained  as  discrete  analogues  of  A,  uqj  and  pQ  j  in  Section  2.4  (sec  [1] 
for  all  the  details). 

3.3.  Multilevel  Algorithms  for  Solving  Problem  (3.11). 

3.3.1.  Synopsis 

We  first  introduce  a  discretization  parameters  h-  to  which  we  associate  all  the  above  discrete 
spaces.  For  simplicity  we  denote  by  l}  the  space  We  assume  that  the  sequence  {ZJ}  satisfies 

the  following  inclusion  property 


(3.12)  Z°  CZ1  c  ...  CZJ. 


At  level  J  (the  finest  level)  we  wish  to  solve  problem  (3.11)  with  h  =  hj. 
v 

Before  defining  a  multilevel  algorithm  for  solving  problem  (3.11),  we  describe  in  the  following 
Section  3.3.2  the  solution  of  general  variational  problems  by  multilevel  methods.  The  application  to 
the  specific  problem  (3.11)  will  be  discussed  in  Section  3.3.3. 

3.3.2.  A  Multi  level  Method  for  Linear  Variational  Problem  in  Hilbert  Spaces. 

Let  V  be  a  Hilbert  space  with  (•  ,  •)  as  inner  product  and  |[-||y  the  corresponding  norm.  We 
consider  the  following  problem 
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Find  ufV  such  that 


(3.13) 

where 


[a(u,  v)=L(v),  VveV, 


(1)  a:  VxV— R  is  bilinear ,  continuous  and  V-elhptic , 

(2)  L:  V— »R  is  linear  and  continuous. 


We  consider  now  a  family  of  finite  dimensional  subspaces  V^CV^C  V“C...C  V^CV.  The 


idea  here  is  to  approximate  (3.13)  by 


(3.14) 


'Find  u^eV^  such  that 
.aj(u^>  v)  =  Lj(v),  V  v€VJ, 


where  aj  and  Lj  are  approximations  to  a(-,-)  and  L  respectively  (for  those  applications  associated  to 
mixed  finite  element  approximations,  aj  and  Lj  are  never  the  restrictions  of  a(-,-)  and  L  to  VxV  and 
V  respectively). 

The  basic  principle  of  all  multilevel  methods  is  to  solve  (3.14)  using  solutions  of  problems  of 
»  j 

the  form  (3.14)  defined  on  V  ,  j  =  0,  1,  ...  J  — 1.  A  classical  way  to  handle  this  is  to  use  a  V-cvcle 
multilevel  method  15,  C,  7,  S].  For  problem  (3.14)  the  V-cycle  with  J  levels  takes  the  following  form: 


Step  0:  Suppose  that  is  known. 

Step  1:  Starting  from  u^,  iterate  i/j  steps  of  some  iterative  method  and  call  the  result  uj^. 

Step  2:  Now  for  j  =  J  —  1,  ....  1,  assuming  that  u  ^  ^  is  known  and  starting  from  0  perform  u  ^  steps 

of  some  iterative  procedure  for  solving  the  following  variational  residual  equation 


15 


(3.15) 


i  J  '  *  J  i 

a-  (uJn,  v)  =  L,(v)-53  »i  Ki  ,  v).  v(V  - 
J  J  1  =  3 

u?fVj. 

V  J 

Call  u*i  the  result  of  this  smoothing. 


Step  3:  For  j=0  solve  exactly  the  residual  equation  (3.15).  Set  u j1  =u®. 

Step  4:  For  j  =  1,  2,  ...,  J,  assuming  uj^  *  is  known,  take  *+unJ  as  an  tniiial  condition 

Perform  fi j  steps  of  some  iterative  procedure  for  solving  (3.15).  Call  the  result  uj^. 


Step  5:  Take  u^,j=u|^. 


3.3.3  •Application  of  the  V-cvcle  Method  to  the  Solution  of  Problem  (3.11). 

Problem  (3.11)  is  a  particular  case  of  problem  (3-14).  Thus,  it  can  be  solved  by  the  multilevel  method 
described  in  Section  3.3.2.  Orce  the  basic  iterative  methods  involved  in  the  V-cvcle  have  been 
specified,  thus  applying  the  above  multilevel  method  is  canonical. 

The  numerical  results  discussed  in  Section  4  have  been  obtained  using  conjugate  gradient  as  a 
smoother  in  Steps  1  and  2,  taking  i/j  =  2.  For  j  =  0  we  also  used  conjugate  gradient  to  obtain  u°.  In 
Step  4  we  employed  one  iteration  of  steepest  descent. 

The  conjugate  gradient  algorithm  for  solving  problem  (3.11)  is  described  in  Section  4  of  [1], 


4.  Numerical  Results 

In  this  section  we  shall  present  the  results  of  numerical  experiments  where  the  mixed 
element/multi-ievel  domain  decomposilion  methods  described  in  Section  2.3  have  been  applied  to  the 
solution  of  test  problems.  The  examples  considered  here  include  both  some  standard  cases  as  well  as 
physical  problems  arising  in  flow  in  porous  media,  such  as  ( 1 . l)-(  1 .3)  of  Section  1.  In  all  our  examples, 
the  discrete  problem  (2.11)  approximating  the  elliptic  problem  (2.1)  has  been  obtained  using  for 
and  V*1  the  Raviart-Thomas  mixed  finite  element  spaces.  A  full  description  of  these  elements  can  be 
found  in  [1]  and  [2] ;  however  for  completeness  we  shall  describe  these  spaces  in  the  following  Section 
4.1. 
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4.1  Mixed  Finite  Element  Appioximations  of  Problem  (2.1). 


Let  ft  be  the  rectangular  domain  (0,  x^)x(0,  y^)  and  let  Ax:  0=Xq<x1<...<xn,  ,=xl  anc^  Ay: 
0=yg<y^<...<y^,  =y^  define  partitions  of  [0,  x^J  and  [0,  VjJ,  respectively.  For  A  a  partition, 
define  the  piecewise  polynomial  space 

Mg(A)  =  {veCs((0,  L]):  v  is  a  polynomial  of  degree  <r  on  each  subirtierval  of  A), 
where  s=  — 1  refers  to  the  discontinuous  functions.  We  define  now  the  following  approximations  of 

f) 

L“(ft),  Il(ft;  div)  and  Vc  respectively 

wJr  =  M |  (Ax.)  ®  Ms  (Ay), 

^'{j  =[_Ms+i(Ax)  ®  M|(  Ay)Jxj^Ms(Ax)  ®  Ms^j(Ay)^J, 

vh,o=vh'n{«  v  ‘'=°  •»  an}, 

where  h  =  max  |(xjj_^—  Xj),  (y._^— y.)j.  remark  that  these  spaces  satisfy 


V-v  r,  V  veV£‘(i.e.  V-V£‘cW£‘) 


,s,r 


s,r. 


rrS,rx 


In  our  numerical  experiments  we  set  r  =  l. 


4.2.  Solution  of  Standard  Test  Problems 

Motivated  by  applications  in  reservoir  engineering  we  are  considering  now  the  following  class  of 
test  problems: 


(4.1)' 


V’(AVP)_f(1:  0)  6(0,  1)> 

AVp-i/  =  Q  or.  5ft, 


where 


ft  =  (0,  1)'  and  where  A 


defined  by  either 
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(0 

or 

(ii) 

or 

(i«) 


A  — A-.  —I, 


a=a2= 


l+100(x2  +  y2) 


I, 


A  =  Ag  =  aI,  where  o  =  100  if  0<x<.5  and  a  =  l  if  .5<x<l. 


The  partitionings  of  Cl  used  to  implement  the  domain  decomposition  are  those  shown  in 
Section  8  of  [1],  In  particular  a  (Nx,  Ny)  decomposition  involves  a  partitioning  into  NxNy  rectangular 
subdomains  whose  edges  are  parallel  to  the  coordinate  axis. 

Table  4.1  depicts  the  number  of  multi-level  V  cycles  versus  mesh  and  subdomain  partitions: 


Coefficient 

h"1 

f— Subdomains.  sS:V  cvcl 

A1 

20 

(4.  6) 

40 

(4,  6);  (16,  7) 

80 

(4,  9);  (16,  8);  (64,  7) 

A2 

20 

(4,  6) 

40 

(4,  8);  (16,  7) 

80 

(4,  10);  (16,  8);  (64,  7) 

'  A3 

20 

(4,7) 

40 

(4,  6);  (16,  7) 

SO 

(4,  10);  (16,  8);  (64,  7) 

Number  of  Cycles  versus  Mesh  Size  and  Subdomain  Partition  for  the  3-Leve!  V-Cycle. 

Table  4.1 


IS 


Interestingly  the  above  table  applies  for  the  three  cases  (i)  — (iii).  We  also  observe  that  the  number  of 
grid  points  by  subdomain  is  the  same  for  the  three  decompositions  considered  and  that  the  number  of 
V  cycles’is  practically  independent  of  h  despite  the  fact  that  the  dimension  of  the  interface  problem  is 
growing  like  h 

To  further  illustrate  the  efficiency  of  the  above  methods  we  are  providing  in  Table  4.2  below 
the  dimensions  of  the  various  finite  element  and  boundary  spaces  involved  in  our  combined  domain 
decomposition/mixed  finite  elements  methodology  (below,  7  is  defined  by  an  N  x  M  decomposition). 


h"1 

Dim  Wh 

Dim  Vh 

Dim  V°, 

oh,  7 

20 

1600 

3120 

40  (N  +  M)  -  79  -NM 

40 

6400 

12640 

SO  (N  +  M)  -159  -NM 

SO 

25600 

50S00 

160  (N  +  M)  -319  -NM 

Dimension  of  the  Discrete  Spaces 
Table  4.2 
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This  insensibility  to  the  smooth  or  fast  variation  of  coefficient  A  over  fi  is  a  remarkable 
property  which  shows  that  this  methodology  has  attractive  potential  for  the  solutio  of  badly 
conditioned  practical  problem  such  as  geostatistics  problems  arising  in  porous  m  -lia  ([9,  10].) 

The  above  results  represent  a  substantial  improvement  in  terms  of  robu- ness  and  speedup 
compared  to  the  results  obtained  in  [1]  for  the  same  test  problems  with  the  same  grids  and 
decompositions. 

Another  interesting  property  of  the  above  methodology  (already  observed  in  [1])  is  that  the 
subdomain  problems  need  not  be  solved  exactly.  We  also  observed,  concerning  the  multilevel  solution 
of  the  matching  problem,  that  one  to  two  V  cycles  are  sufficient  in  practice  to  achieve  the  solution 
within  truncation  error;  in  particular,  with  u^  =  fj^  =  2  in  the  algorithm  of  Section  3.3.2,  the  initial 
residual  is  reduced  by  six  orders  of  magnitude  in  six  to  seven  iterations,  the  largest  reduction  taking 
place  in  the  first  V-cycle. 

4.3.  Solution  of  Real-Life  Test  Problems. 

To  be  honest  the  test  cases  discussed  here  are  more  relevant  to  [1]  since  the  domain 
decomposition  methodology  is  exactly  the  one  described  in  the  above  reference,  i.e.  without-yet- 
multilevel  speedup.  Nevertheless,  we  have  inserted  these  problems  because  they  are  typical  of  real-life 
applications  in  petroleum  reservoir  engineering.  Also  they  provide  significant  benchmarks  for  elliptic 
solvers  of  various  types. 

This  first  problem  to  be  considered  was  communicated  to  us  by  petroleum  reservoir  engineers. 
It  is  a  model  for  a  discrete  shale  barrier  and  involves  solving  (1.1)-(1.3)  where  A  is  visualized  in  Figure 
4.1,  where  we  have  used  different  scales  for  L  and  H  since  L  is  of  the  order  of  300  feet  and  H  is  of  the 
order  of  20  feet  implying  an  aspect  ratio  of  15.  Also  the  thickness  of  the  barrier  is  of  order  one  foot. 
The  ratio  of  permeability  coefficients  is  10“. 
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DISCRETE  SHALE  BARRIER  PROBLEM 


v 

Geometry  of  the  Discrete  Shale  Barrier  Problem 
Figure  4.1 

The  arrows  in  Figure  4.1  indicate  the  flow  direction. 

.  Concerning  the  numerical  solution  of  the  above  problem  we  have  been  using  a  40x40  finite 
element  grid  and  a  (2,  2)  domain  decomposition.  For  comparison  purposes  we  have  treated  the  cases 
with  aspect  ratios  1  and  15. 

Using  the  domain  decomposition  algorithm  discussed  in  [1]  we  need  33  iterations  if  R=1  and 
48  if  R=15.  We  can  expect  the  number  of  iterations  to  be  practically  independent  of  R  once  our  V 
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In  the  same  vein  the  second  problem  is  also  a  real  life  problem  (l.l)-(l .3)  where 
A=k(x,  y)/p(c)  and 

/j(c)=cMr1/4+(i~c)Mr1/4> 

with  pj,po>0. 

Applying  the  domain  decomposition-mixed  finite  element  methods  of  [1]  to  the  above  problem, 
with  a  80  x  80  finite  clement  grid  and  a  (10,  10)  domain  decomposition,  the  solution  was  obtained  in 
9  conjugate  gradient  iterations.  This  represents  a  substantial  improvement  over  a  preconditioned 
conjugate  gradient  solution  of  the  same  discrete  problem  (without  domain  decomposition)  since  the 
convergence  was  requiring  then  about  150  iteration,  (taking  advantage  of  a  good  unital  guess). 
Incidentally  the  lowest  order  Raviart-Thomas  space  (r  =  0  (in  4.1))  or  cell-centered  finite  differences  [11] 
do  not  work  well  on  this  type  of  problems  due  to  the  impossibility  for  these  low  order  approximations 
to  reproduce  correctly  flows  which  are  not  parallel  to  the  coordinate  axes;  this  drawback  disappears  if 
we  chose  r=l. 

In  Figure  4.2  we  have  visu  .»ized  the  permeability  k(x,  y),  this  data  was  measured  by  researchers  at 
Atlantic  Richfield  Corporation  and  kindly  communicated  to  us.  Similarly  the  function  A=k/p  is 
visualized  in  Figure  4.3. 

5.  Mesh  Refinements  Via  Domain  Decomposition 

Mesh  refinements  are  necessary  when  strong  gradients  arise  locally.  In  view  of  saving  computer 
storage  and  avoiding  complicated  data  structures  it  is  interesting  to  incorporate  local  grid  refinement 
over  subdomains  where  the  strong  variations  are  arising  and  retain  coarser  grids  elsewhere.  The 
concept  of  domain  decomposition  provides  an  elegant  and  systematic  way  to  implement  the  above 
ideas.  In  this  section  we  would  like  to  present  a  particular  implementation  of  our  scheme,  new  to  our 
knowledge,  relying  again  on  a  combination  of  Raviart-Thomas  mixed  finite  element  and  domain 
decomposition  methods. 


'  '  '•  •  *’  •  '■  ;‘\  '  *  ■  i 

f  ■  ;  ■  •  -  -  -  -  ‘  •  - 


i*' 


Representation  of  k(x,  y) 
Figure  4.2 


Representation  of  A  =  k//r. 
Figure  4.3 
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5.1  Mesh  Refinement  Via  a  Modified  Raneart -Thomas  Mixed  Finite  Element  Method 


Consider  the  situation  depicted  in  Figure  5.1  where  a  local  refinement  is  necessary  m  a 
subregion  fi*  of  Q.  The  basic  idea  is  to  employ  essentially  mixed  finite  elements  of  Raviart-Thomas 
type  inside  and  outside  subregion  0*;  the  main  issue  here  is  clearly  the  matching  between  the  ’’fine" 
and  ’’coarse”  approximations.  To  realize  this  matching  we  introduce  the  following  finite  dimensional 
spaces  of  mixed  type. 

Let  Q  =(a  ,  b  )x(c  ,  d  )  and  define  A£  and  Ay  be  partitions  of  [a*,  b*]  and  [c*,  d*], 
respectively.  Generalizing  the  notation  of  Section  4.1,  we  denote  by 


(5.1) 

W1  r* 

V 

(5.2) 

v*-1- r 

h* 

and 

(5.3) 

V*-1.  * 

h*.0 

(ft*)  =  Mf*(Ax*)  ®  Mr*(Ay*), 

(G*)  =  (Mj  +  1(Ax*)  ®  M.r*(Ay*))  x  (M^j (Ax*)  ®  mJ 

(Q  )  =  V^  ’  (fl*)  D  {q  :  q-i/  =  0  on  3fl*}. 


(Ay*)), 


Similarly  we  define  the  corresponding  "coarse”  spaces  by 

(5.4)  w^1,  r(n-o*)=wj;’  r| 

n-n* 
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Figure  5.1 


and 


(5.5) 

v;1,  r(n-n*)=v;1,  r  i 
h  ;  i»  'n-n*’ 

with  W^’ 

r  and  vj’  r  as  defined  in  Section  4.1.  We  set 

(5.6)  : 

wh  =wh* r  (-Q*)  u  w;1’  r(n-n*), 

(5.7) 

V 

vh  =vh*  r*(Q*) u  vh’  r(n-«*). 

(5.8) 

Vh,  0  =  Vh  n  S-^  =  0  on  5f2}. 

Strictly  speaking  \\f  and  vf  are  not  Raviart-Thomas  spaces,  however,  they  share  the  same 

approximating  properties  which  include  div  Cwf  and  the  order  approximation  is  the  same  if 
r*=r. 

.  From  a  computational  point  of  view  this  refinement  technique  is  well  suited  for  domain 
decomposition  with  Q*  and  as  subdomains. 
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The  above  approach  is  well  suited  for  a  multi-level  solution  of  problem  (2.1)  in  which  we  shall 
use  different  number  of  grid  levels  in  the  subdomains  (usually  more  grid  levels  in  the  more  refined 
regions).'  Domain  decomposition  allow  a  lot  of  flexibility  by  the  fact  that  in  one  of  the  phases  of  their 
realization  they  decouple  the  computation  to  be  done  in  each  subdomain. 

6.  Conclusions 

From  the  numerical  results  described  in  this  paper  the  combination  of  mixed  finite  element, 
domain  decomposition  and  multilevel  methods  discussed  in  Sections  2,  3  and  4  provides  a  robust, 
accurate  and  fast  technique  for  solving  elliptic  problem  with  non-smooth  coefficients  like  those  arising 
in  flow  in  porous  media  and  other  applications  from  Mechanics  and  Physics. 

These  methods  are  quite  interesting  from  a  parallel  computing  point  of  view  since  the  ratio 

Work  in  Solving  Subdomain  Problems 
Communication  Costs 


is  of  order  0(h  *). 

Here  the  communication  involves  the  transfer  of  the  boundary  data  at  the  subdomain  -terfaces. 

We  are  presently  cooperating  with  the  computer  scientists  at  the  National  Science  Foundation 
Center  for  Research  in  Parallel  Computation  in  the  parallel  implementation  of  the  methods  discussed 
in  this  paper. 
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ABSTRACT 

The  main  goal  of  this  paper  is  to  discuss  mixed  variational  formulations  for 
time  dependent  problems  such  as  initial  and  boundary  value  problems  for  the  heat 
and  wave  equations  in  a  bounded  domain  fi  of  H^(JV  >  l).  Then  we  shall 
use  these  formulations  to  derive  mixed  finite  element  approximations  of  the  heat 
and  wave  equations.  Finally,  an  application  to  an  exart  boundary  controllability 
problem  for  the  wave  equation  will  be  presented  together  with  some  numerical 
results.  The  techniques  and  application  briefly  considered  here  will  be  discussed 
with  more  details  in  a  forthcoming  paper. 

INTRODUCTION 

Mixed  variational  principles  and  the  associated  finite  element  approximations 
have  proved  to  be  very  useful  in  order  to  derive  accurate  solution  methods  for 
boundary  value  problems  for  partial  differential  equations.  This  is  particularly  true 
for  elliptic  problems  (see,  e.g.,  [l],  [2]  and  the  references  therein)..  A  strong  point 
of  these  techniques  -  compared  to  more  traditional  finite  element  methods  -  is  that 
they  give  fairly  accurate  approximations  of  the  derivatives ;  this  iast  property  is 
very  interesting  since  in  many  problems  one  is  more  interested  by  the  derivatives 
of  a  function  than  by  the  function  itself.  Mixed  methods  have  also  been  applied 
to  time  dependent  problems  (see,  e.g.,  { 3] )  but  there  are  indeed  very  few  published 
papers  and  applications  where  these  methods  have  been  used  for  time  dependent 
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problems  compared  to  the  more  classical  finite  element  methods.  Motivated  by 
optimal  control  applications  (cf.  [4],  [5])  we  shall  briefly  discuss  in  this  short  article 
the  following  topics: 

(i)  Mixed  variational  formulations  for  the  heat  and  wave  equations  (Section  1.). 

(i>^  Mixed  finite  element  approximations  of  the  heat  and  wave  equations  (Section 

2.). 

(iii)  An  application  to  a  boundary  control  problem  for  the  wave  equation  (Section 
3.). 

1.  MIXED  VARIATIONAL  FORMULATIONS  FOR  THE  HEAT  AND  WAVE  EQUATIONS. 

1.1  Formulation  of  the  basic  time  dependent  problems. 

Let  fl  be  a  bounded  domain  of  JEtN(N  >  1)  ;  we  denote  by  T  the  boundary  of 
D  .  Let  T  be  a  positive  number  (possibly  equal  to  +oo)  ;  we  denote  by  Q  and  E 
the  following  sets  of  R  f+1  : 

Q  =  0  x  (0,T),E  =  T  x  (0,T). 


We  suppose  now  that  physical  phenomena  are  taking  place  on  fl  ,  modelled  by 
either  the  following  heat  equation 


(1.1) 

u*  -  Au  =  /  t'n  Q, 

(1.2) 

u  =  g  on  E, 

CO 

u(x,0)  =  u0(x)  on 

or  by  the  following  wave  equation 


vtt  -  Au  =  /  in  Q, 
u  =  g  on  E, 


u(x,  0)  =  u„(x),  ut(x,  0)  =  uj  (r)  on  fl. 


(1.4) 

(1.5) 

(1.6) 


O 


In  (l.l)  -  (1.6)  we  have 


/  \ff  du 

*  =  {*<>,"  i,«*  =  apu«  = 


—  A  _  Vil 

at2’  f-ra  if’ 

t  =  l  1 


It  follows  from,  e.g.  [6],  [7],  that  each  of  the  two  above  problems  has  a  unique 
solution  provided  that  the  data  /  and  $  belong  to  well  chosen  functional  spaces. 
Since  this  paper  is  engineering  oriented  we  shall  not  go  into  the  details  of  those 
(Sobolev  type)  spaces  for  which  the  above  problems  are  well-posed  (there  will  be 
however  some  exceptions). 

1.2  Mixed  variational  formulations  for  problems  (1  )  —  (1.3)  and  (1.4)  —  (1.6). 

The  key  idea  is  to  take  Vu(V  =  {af-}£Li)  as  master  variable  ;  we  introduce 
therefore  a  new  unknown  p  defined  by 


(1.7) 


p  =  Vu(m  Q). 


Assuming  that  u  and  p  are  sufficiently  smooth  we  obtain  -  integrating  by  parts  with 
respect  to  the  space  variables  -  the  following  mixed  variational  formulations: 
Mixed  variational  formulations  of  the  heat  equation  (l.l)  —  (1.3)  : 

(1.8)  f  (ut  -  V  -  p  -  f)vdx  =  0,  VveL2(fi),a.e.  on  (0,T), 

J  o 

(1.9)  j  (p  •  q  +  uV  •  q)dx  =  I  gq  •  ndT,  VqeH (H,  div),  a.e.  on  (0,  T), 

J  n  Jr 

(1.10)  u(x,  0)  =  u0(x)  on  H. 

Mixed  variational  formulations  of  the  wave  equation  (l.-i)  —  (1.6)  : 

(1.11)  j  (utt  —  V  •  p  —  f)vdx  =  0,  VueL2(n),  a.e. on  (0,T), 

Jn 

1.12)  J(  p.  q  +  uV  •  q )dx  =  J  gq  •  ncT,  VqehT(Q,  div),  a.e.  on  (0,  T), 

(1.13)  u(x,  0)  =  U0(x),Ut( 2,0)  =  Ui(2).D 

In  (1.8)  -  (1.13),  we  have  used  the  following  notation:  y  •  z  =  Vy, zelR^jn 

is  the  unit  vector  of  the  outward  normal  at  T;dx  =  dxi  ■  ■  ■  dx //  and  finally 

H{ Cl,  div)  =  (q|qeL2(n),  V  ■  qfL2(n)}. 


3 


2.  MIXED  FINITE  ELEMENT  APPROXIMATIONS  OF  THE  HEAT  AND  WAVE  EQUATIONS. 
2.1  Generalities. 

With  h  a  space  discretization  step,  we  approximate  L2(fl)  and  H(Cl,div)  by  Vh 
and  Qh.  ,  respectively.  We  suppose  that  Vh  C  L7{Cl),Qh  C  H( U,  div)  and  also  that 
Vh  and  Qh  satisfy  compatibility  conditions  implying  convergence  properties  for  the 
corresponding  approximations  (see  e.g.,  [l],  (2]  for  details);  an  important  condition 
to  be  satisfied  is: 

(2-1)  V  •  Qk  C  VK. 

In  the  particular  case  where  fl  is  a  2  dimensional  polygonal  whose  boundary  is  the 
union  of  segments  parallel  to  the  coordinate  axis,  we  associate  to  fl  a  “partition” 

Rh  such  that 

(•')  A  =  W.n  =  KuiR~t ?, 

(ii)  Each  K  is  a  rectangle  whose  edges  are  parallel  to  the  coordinate  axis, 

(iii)  If  K  and  K'eRh ,  then  K  C\  K'  =  <p,  and  either  K  n  K '  —  <f>,  or  K  and  K '  have 
only  a  whole  edge  or  one  vertex  in  common. 

Following  [l],  [2]  and  [8]  -  (lOj,  a  convergent  choice  for  Vh  and  Qh,  constructed 
from  the  above  Rh,  is  given  by: 

(2-2)  Vh  =  {vh\vh\KtQkyKeRh}, 

f  Qh  =  {qa-q*  =  {q</i}?=i.q*i*e(-P*+i  ®  Rk)  x  (p*  ®  P*+i), 

(2.3)  <  VKeRhiqih  is  continuous  along  the  edges 

l  of  Rh  parallel  to  0£,  + i},' 

in  (2.2),  (2.3),  k  is  a  nonnegative  integer,  Qk  =  Pt  0  P*,P,  is  the  space  of  the 
polynomials  in  one  variable  of  degree  <  s,  and  :  +  1  has  to  be  taken  modulo  2. 

With  such  a  choice  for  Vh  and  Qh  ,  condition  (2.1)  is  clearly  satisfied. 
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2.2  Discretization  of  the  heat  equation  (1.1)  -  (1.3). 

Semi  —  Discretization  in  space  : 

Using  the  spaces  Vj*  and  Qh  we  shall  “space  discretize”  (1.1)  -  (1.3),  via  (1.8) 
-  (1.10)  as  follows: 

Find  a  pair  (ufc(t),Pfc(t)}cV^  x  Q a.e.  on(0,  T),  such  that 

(2-4)  ^(^T  ~  V  ‘Ph  ~  fdVhdx  -  °yvh£Vh,  a.e.  on  (0 ,T), 

(2.5)  [  (p*  •  qa  +  u/,V  •  qjdx  =  f  ghqh  •  ndT;VqheQh,  a.e.  on  (0,T), 

Jo  Jr 

(2.6)  Ufc(0)  =  u0*. 


In  (2.4)  -  (2.6),  fh,gh  and  u0h  are  convenient  approximations  of  f,g  and  uQ,  respec¬ 
tively  (we  can  take,  for  example,  uQh  as  the  jL2-projection  of  u0  on  Vh)  . 

The  above  approximation  is  not  practical  since  we  still  have  to  solve  an  ordinary 
differential  system,  or  to  be  more  precise  a  system,  coupling  ordinary  differential 
equations  and  (linear)  algebraic  equations. 

Full  Discretization  in  space  —  time  :  Concentrating  (for  simplicity)  on  the  back¬ 
ward  Euler  time  discretization  of  (2.4)  -  (2.6)  we  finally  obtain  the  following  system 
of  difference  -  algebraic  equations  (with  Af(>  0)  a  time  discretization  step): 

For  n  >  0,  find  {u£+1,  p£+1}eVk  x  Qh  such  that 


(2.7) 

(2.8) 

(2.9) 


ttfc  =  “oh, 


u 


n-rl 


—  U. 


'n 


At 


-  v  •  Pr 1  -  vhdx  =  0,VvhcVH,  • 


/  (p£+1  •  qh  ~  ■  qh)  dx  =  [  ?£+1qA  •  ndr.Vq^cQ*. 

Jo  Jr 


From  a  practical  point  of  view,  we  can  easily  eliminate  u£+1  from  (2.8),  using 
the  fact  that  V  •  q^tV/j  ;  we  obtain  then  the  following  linear  variational  equation 
satisfied  by  p£+1  : 

/n(AtV  •  Pr*V  •  q^  +  p£+1  ■  q  h)dx  =  Jr  g£+1  qh  ■  nd T 
-  /n(u-  +  At/;  +  1)V  •  qhdx,\fqheQh]p^hQh. 


(2.10) 
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Solving  (2.10)  can  be  done  by  a  direct  method  -  such  as  Cholesky’s  since  the  bilin¬ 
ear  form  in  (2.10)  is  symmetric  and  positive  definite  -  or  by  a  conjugate  gradient 
algorithm  (see,  for  example,  [ll]).  Once  p£+1  is  known,  computing  u£+1  from  (2.8) 
is  straightforward. 

Similarly,  instead  of  backward  Euler,  we  could  have  used  schemes  such  as  for¬ 
ward  Euler,  Crank  -  Nicholson,  multisteps,  Runge  -  Kutta,  .... 

2.3  Discretization  of  the  wave  equation  (1.4)  —  (1.6). 

Starting  from  the  following  variant  of  (2.4)  -  (2.6):  Find  a  pair 


{u/l(t),pA(f)}cV1il  x  Qh,  a.e.on(0,T),  such  that 


f  <9^ 

(2.11)  fh)vhdx  =  0,VvheVh,  a.e.  on(Q,T), 

(2.12)  /  (Ph  •  q/i  +  uh V  •  qK)dx  =  f  ghqh  *  ndT,  Vq^eQ*,  a.e.  on(0, T), 

Jn  Jr 

(2.13)  uA(0)  =  uo„^(0)  =  u1A, 


we  can  fully  discretize  the  wave  problem  (1.4)  -  (1.6)  by  the  following  variant  of  the 
usual  second  order  accurate,  explicit  finite  difference  discretization  scheme  of  the 
wave  equation: 

Assuming  that,  for  n  >  0,  u£,p£  and  u£-1  are  known  compute  first  <+>  as 
the  solution  of 


(2.14) 


-  V-p £  -  fff)vhdx  =  0 ,VvheVh-,u£+1eVh, 


and  then  p£+1  as  the  solution  of 

(2.15)  [  p^+1-q,Hdx=  [  gZ+1qh-ndT-  [  u£+1  V  •  qhdx,  Vq^cQ,,; p£+1eQh. 

Jn  Jr  Jn 

A  most  important  step  is  clearly  the  initialization  of  scheme  (2.14),  (2.15);  assum¬ 
ing  that  f,g,u0,u i  are  sufficiently  smooth  we  shall  proceed  as  follows:  compute 

and  uh 


(2.16) 


U/t  —  u0h,  u^  —  uh  +  2Afui/j, 
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(2.17) 


p  fcQh, 

/n  P h  *  <lhdx  =  Jr  •  ndT  -  /n  uohV  •  qAdx,  Vq^eQ*. 


As  shown  in  (12],  Uh(t)  and  Ph(t)  will  converge  to  u(t)  and  Vu(t)(u  :  solution  of 
(1.4)  -  (1.6))  as  h  and  At  — ►  0  if  a  stability  condition  such  as 

(2.18)  At  <  Ch 

is  satisfied. 

Second  order,  unconditionally  stable  implicit  variants  of  the  above  scheme  can 
be  obtained;  they  will  discussed  in  a  following  paper,  together  with  applications  to 
boundary  control  of  the  wave  equation. 

3.  APPLICATION  TO  AN  EXACT  CONTROLLABILITY  PROBLEM  FOR  THE  WAVE 
EQUATION,  VIA  DIRICHLET  BOUNDARY  CONTROLS. 

3.1  Formulation  of  the  boundary  control  problem. 

We  follow  here  [4],  [5];  we  consider  then  a  phenomenon  taking  place  in  Cl  and 
modelled  by  the  wave  equation  (we  keep  the  notation  of  Section  1): 

(3.1)  utt  —  A u  =  0  in  Q , 
with  the  initial  conditions 

(3.2)  u(x, 0)  =  u0(x),ut[x,o)  =  ui(x)  in  Cl. 

The  problem  here  is  to  find  g  defined  over  ^  x  (Oj^))  such  that  the  following 
final  conditions 

(3.3)  u(x,T)  =  0,  ut(x,T)  =  0  on  Cl 
hold  if  one  has 


(3.4) 


u  =  g  onE 
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as  boundary  condition. 

It  has  been  proved  by  several  authors  (see  [4j,  [5],  [l3j  for  references)  that  such  a 
g  exists  provided  that  T  is  sufficiently  large  (the  lower  bound  of  the  7”s  for  which 
(3.3)  holds,  Vuo.ui,  is  -  not  surprisingly  -  of  the  order  of  diameter  {Cl)). 

3.2  Calculation  of  an  exact  Dirichlet  control  via  the  HilbcrtUniqucncssMcthod 
of  J .  L.  Lions 

In  [4],  [5],  J.L.  Lions  has  introduced  and  analyzed  a  systematic  way  for  con¬ 
structing  Dirichlet  controls  for  which  (3.3)  holds.  The  construction  technique  is 
systematic  and  based  on  the  Hilbert  Uniqueness  Method  (HUM)  to  be  briefly  dis¬ 
cussed  below.  From  now  on,  we  suppose  that 

(3.5)  tI„ei2(f!),ul£^-1(n)(=  (ff0‘(n)'), 

where 

/j 

Hl0(Cl)  =  {v\vtL*{n),  -~eL2{Cl)t\/i  =  1, •  •  •  N, v  =  0  on  T}, 

ox% 

H~l{Cl)  is  the  dual  space  of  H^{Cl), 
and  we  define  E  and  E'  by 

(3.6)  e  =  Hl{ n)  x  l2{cl),e'  =  tf-^n)  x  L2(n). 

Next  we  define  an  operator  AeL(E,  E')  as  follows: 


(0 

Take  e  =  {e0,  ei}e£; 

(«•*■) 

Integrate  from  0  to  T  : 

(37), 

<f>tt  —  Ad>  =  0  inQ, 

(3.7), 

4>  =  0  on  y~], 

(3.7)3 

d>(x,0)  =  eo(x),^i(z,0) 

(m) 

Integrate  from  T  to  0  : 

(3.8), 

xfrtt  —  At/>  =  0  inQ , 
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(3.8)2 

dif> 

(3.8)3 

i{>(xtT)  =  Q,rl>t[x,T)  =  0  onfi. 

(»u)  take 

(3.9) 

Ae  =  (*(0),-*(0)}, 

where  tp(0)  ( resp .  V’t(O)  )  is  the  function  x  —*  ip(x,  0)  (resp.  x  —*  ipt(x,0)). 

It  follows  from  J.L.  Lions  [4],  [5]  that  A eL{E,E')yT  >  0;  moreover,  if  T  is 
sufficiently  large  ( T  >  diameter  (fl)  )  then  A  is  a  strongly  elliptic  operator  from  E 
onto  E\  In  addition  to  these  properties,  A  is  self-adjoint  and  satisfies  (with  obvious 
notation): 

(3.10)  (A e,e/)=/*  ^^dTdt, Ve,e'e£; 

E 

in  (3.10),  (*,•)  denotes  the  duality  pairing  between  E'  and  E  which  satisfies 

(Ae,e')  =  f  (Ae)-e'dx 
J  n 

if  Ae  is  sufficiently  smooth. 

Application  to  the  exact  boundary  controllability  of  the  wave  equation  : 

(i)  Solve 

(3.11)  Ae  = 

(ii)  Solve  (3.7),  taking  for  e  ,  in  (3.7)3  ,  the  solution  of  (3.11). 

(iii)  Take  g  =  on  £ . 

If  T  is  sufficiently  large,  it  follows  -  from  the  properties  of  A  -  that  (3.11)  has  a 
unique  solution  in  E  ;  we  have  (cf.  [4],  [5])  gtL2{Yff),  and  the  corresponding  solution 
of  (3.8)  satisfies  (3.1)  -  (3.4),  implying  that  g  is  a  Dirichlet  boundary  control  for 
which  the  exact  controllability  property  (3.3)  holds.  Actually,  of  all  the  Dirichlet 
boundary  control  for  which  exact  controllability  holds,  the  one  obtained  by  HUM, 


9 


i.e.  by  solving  (3.11)  is  the  only  one  of  minima/  norm  in  Z-2(5I),  as  shown  in  [4], 
(5j.  From  the  properties  of  A  ,  problem  (3.11)  can  be  solved  by  a  conjugate  gradient 
algorithm  operating  in  space  E  ;  such  an  algorithm  is  described  in  [l3j,  [14],  together 
with  conforming  finite  finite  element  implementations  of  it. 


3.3  Mixed  formulation  of  the  boundary  control  problem. 

In  fact,  we  shall  describe  a  mixed  formulation  of  problem  (3.11): 

Assuming  that  the  initial  data  ua  and  ui  are  sufficiently  smooth,  so  that  we 
can  use  integral  representations,  the  problem  is  now  to  find  a  triple  {e0,p0,  ej) 


satisfying 

(  {e0,Po}£^o,eieL2(n);V{u0,7r0}£W<>,t;1eL2(n)  we  have 

[Ot  1.2  j  \ 

l  Jn&ti °)v<>  ~  V>(0)vi )dx  =  Jn(u !Va  -  u0vi)dx , 
where  in  (3.12): 

(i)  The  space  W0  is  defined  by 

(3  13)  /  W°  =  {(l,o,ff0}|woe^2(n),7r0c(L2(n))wr  fn(7r0-q  + v0V-q)dx  =  0, 

I  Vqeii(n,d:t;)}; 

it  can  be  shown  that 


{  V0  f  ^ O  }eW0  ~  VotHl{n),*o  =  Vu0. 

(ii)  t/>(0)  and  V’t(O)  are  obtained  from  e0,p0,  e\  as  follows: 

Integrate  from  0  to  T  the  mixed  formulated  following  wave  equation  (cf.  Section 


2): 

(3.14)! 

/  (^tt  —  \  •  p)udx  =  0,VueL2(n),  a.e.  on{0,T), 

J?t 

(3.14)o 

/  (p  •  z  —  <?V  •  z)dx  =  0,  Vz£i7(f2, div),  a.e.  on  (0,T), 
Jn 

(3.14)3 

<f>{x,0)  =  s0iz),<pt(x, 0)  =  ei(x)  on  f7; 

then  from  T  to  0  (using  the  fact  that  ~  —  p  •  n  on  Tjj): 

(3.15)  x  f  (vti  —  V  •  qK'fiz  =  0,VveL2(n),  a.e.  on(0,T), 

Jn 

(3.15) o  /  (q  ■  z  4-  •  7.)ci  =  j  p  •  n  z  •  ndr,  Vze//(n,  div),  a.e.  on(0,T), 

Jn  Jr 

(3.15) 3 7")  =  =  0  on  ft. 
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An  easy  calculation  will  show  that  (with  obvious  notation): 

(  fn(M°)eo  ~  H°)ei)dx  =  I  P  *  nP'  ’  ndT dt, 

(3.16)  <  E 

l  V{e0,7r0;ei},{e;,7r '0\e\}tW0  x  L2(Cl). 

From  (3.16)  it  appears  that  the  bilinear  form  occuring  in  (3.12)  is  symmetric  and 
positive  semi  definite  ;  actually,  for  T  sufficiently  large  it  is  strongly  elliptic  (co¬ 
ercive)  over  (W0  x  L2(fi))2.  From  these  properties,  problem  (3.12)  can  be  solved 
by  a  conjugate  gradient  algorithm  operating  in  Wa  x  L2(Cl)  ;  such  an  algorithm  is 
described  in  Section  3.4. 

3.4  Conjugate  gradient  solution  of  problem  ( 8.12 ). 

3.4.1.  Generalities. 

Problem  (3.12)  is  a  particular  case  of 


(3.17)  Find  ueV  such  thata(u,v)  =  L(v),VveV, 


where  in  (3.17): 

(i)  V  is  an  Hilbert  space ,  equipped  with  the  scalar  product  (-,  ■)  ,  and  the  corre¬ 
sponding  norm  1 1  •  1 1  . 

(ii)  a  :  V  x  V  — »  R  is  bilinear,  continuous  and  V-  elliptic  (i.e.  3a  >  0  such  that 
a(«,v)  >  a||t;||2, VueP). 

(iii)  L  :  V  — ♦  R  is  linear  and  continuous. 

It  is  well  known  (cf.,  e.g.,  [15,  Appendix  1])  that  under  the  above  hypotheses, 
problem  (3.17)  has  a  unique  solution.  If  in  addition  to  (i)  -  (iii),  the  bilinear  form 
a(-,  •)  is  symmetric  then  problem  (3.17)  is  equivalent  to  the  following  minimization 
one 


(3.18) 


ucV, 

J(u)  <  J(v),VveV, 


with  J(v )  =  |a(u,u)  —  L{v).  Problem  (3.17),  (3.18)  can  then  be  solved  by  the 


following  conjugate  gradient  algorithm: 


11 


Initialization 


(3.19) 


u°eV  is  given. 


Solve  then 


(3.20) 


(9°eV, 

l  =  a(u°,v)  -  L(v), 


VveV. 


If  9°  —  0,  or  is  “ small ",  take  u  —  u°\if  not,  set 


(3.21) 


w°  = 


Now  for  n  >  0  ,  suppose  that  un,gn,wn,  are  known  with  wn  0  ;  define  then 
un+1,gn+1,wn+1  as  follows: 

Descent :  Compute 


(3.22) 


Pn  =  |kn||2M^n^n), 


(3.23) 


un+1  =  un  —  pntun. 


Test  of  the  convergence  and  updating  the  descent  direction  :  Solve 


{[9  \v)  =  (gn,v)  -  pna(wn,v),\ 

U  <7n+1  =  0  -  or  is  small-take  u  =  un+1  ;  if  not  compute 


(3.25) 


7n  =  |!?n+1i|2/|lffn||2» 


and  update  wn  by 


(3.26) 


wn+1  =  gn+1  +~inwnn 


Do  n  =  n  +  1  and  go  to  (3.22). 


The  above  algorithm  converges,  Vu'W,  and  we  have  (cf.  [16]): 


(3.27)  ||«*-«t|<C|K-u||(^fy)*, 
where  C  is  a  constant,  and  where  the  condition  number  ua  is  given  by 

(3.28)  Ua~  l*sa(V'V^veSa(V'V^' 

with  S  =  {v\veV,  ||t/||  =  l}. 

3.4.2 

Since  problem  (3.12)  is  a  particular  problem  (3.17),  with  V  =  WQ  x  £,2(fl),  it 
can  be  solved  by  the  conjugate  gradient  algorithm  (3.19)  -  (3.26).  An  important 
practical  issue  is  the  proper  choice  of  the  scalar  product  to  be  used  over  WQ  x  L2(fi). 
A  fairly  convenient  one  is  provided  by 


Application  to  the  solution  of  problemjS  .12) 


(3.29) 


fn(v°vo  +  7ro‘K  +  vi  v[)dx, 

v{v0,7r0;vi},{vo>K'>vi}elvo  x  £2(n). 


Applying  algorithm  (3.19)  -  (3.26)  to  the  solution  of  problem  (3.12),  with  W0xL2(Cl ) 
equipped  with  the  scalar  product  (3.29),  we  obtain  the  following  algorithm: 
Initialization  : 


(3.30)  {e°,p°}t\V0,e°eL2(n)  are  given. 

Integrate  then  from  0  to  T  the  wave  equation 


(3.31)  i 

(3.31) 2 

(3.31) 3 


[  [o?t  —  V  •  p°)vdx  =  0,  VveL2(Cl),  a.e.  on  (0,  T), 

J  o’ 

j  (p°  •  z  +  <fCjV  •  z)dx  =  0, Vze/f(fl, dtv),  a.e.  on  (0,T), 
Ai 


<?°(o)  =C^°(o) 


e 


o 
1  • 
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Then  from  T  to  0  : 


(3.32)  i  [  (V&  -  V-q>dx  =  0,Vv£L2(n),  a.t.  on  (0,T), 

J  o 

3.32)2  /  (q 0  -  z  +  rp°V  '  z)dx  =  [  p°  •  n  z  •  ndr,Vz£tf(n,dit;), 

Jn  Jr 

a.e.  on  (0, T), 

(3.32) 3  t/,°(T)=0,tf(T)=0. 

Compute  then  {<7®,  tt<7£}  and  5®  as  follows: 

Solve  the  mixed  elliptic  problem: 

Find  {g%,Trg%}zW0  such  that 

(3.33)  j  [  {g°0  -  V  •  7rff>dx  =  f  (V’f(O)  -  tnjvdx, Vt;eL2(n), 

Jn  Jn 

(3.33) 2  [  (vrgS  ■  q  +  g£V  ■  q)dx  =  0,VqeH(n,div), 

J  n 


and  then 


(3.34) 


0?  =  ti0  -  ^°(0). 


//  {ffo,7r<7°}  =  {0,0}, g£  =  0  ,  or  are  small,  take  p°  •  n|  as  boundary  control;  if 

u 

noi,  set 

(3.35)  {w°0,  zw°\  tu?}  =  {g°,  ng°]  g?}.D 

Then  for  n  >  0  ,  assuming  that  p£),  e?,  <^n,  0n,  {g£,  7rg£},  g?,  {u;£,  ttu>£},  tr? 
arc  known,  we  compute  {e£T^  p£+1},  e?+\ <£n+\ i/>n+1,  {g£+1, 7rg£+1,  }, 

+11  {u-’o"!’1>  -U>£+1},  tu”+1,  <35  follows  : 

Descent  : 

Inteyiute  from  0  to  T 

(3.36) !  f  (<?tn£  —  V  •  pn)vdx  =  0,  Vt;eZ,2(n),  a.e.  on  (0,  T), 

Jn 

(3.36) 2  [  (p’1  •  z  +  •  z)dx  =  0,  VzeJ?(n,  div),  a.e.  on  (0 ,T), 

J  n 

(3.36) 3  ^(0)  =  ^,^(0)  =  <. 
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Then  from  T  to  0  : 


f  {fft  -  V  •  q n)vdx  =  0,Vu£L2(n), 

J  o 


(3.37) i  a.«.  on  (0,T), 

j  (If1  •  z  +  •  z)dz  =  f  p*  •  nz  •  ndl\  Vzc.ff  (n,d:v), 

./n  ir 

(3.37) 2  a.e.  on  (0,T), 

(3.37) 3  f(r)  =  o,v?(r)  =  0 

So/t>e  nou;  the  mixed  elliptic  problem  : 

Find  7Tffo}£^/o  sucA  *Aai 

(3.38)  x  [  (g%  -  V  •  Ttg2)vdx  =  /  V?(0)t;dz,  VveL2(n), 

J  n  do 

(3.38) 2  [  (ng2  ■  q  +  g^V  ■  q)dx  =  oyqeH{n,div), 

in 


and  set 


(3.39) 


Compute  now 


(3.40) 


and  tAen 


(3.41) 

(3.42) 

(3.43) 
(3-44) 


T(o). 


=  L(|fl?|a  +  l»rg?|8  +  lgri,)«to 

(0)-tojkv'*(0))<i* 

jniyou-?-r*eZ-*w?+v?w?)dx 


{<+\pS+1.*?+:}  =  (Cp2.«i 
{^n+i,pa+i}  =  {<?n,pn}-Pn{r,r}, 
{^n-,‘1,qn+1}  =  {v-n,qn}  -  Pn{r,q“}, 
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Test  of  the  convergence.  New  descent  Direction  : 

If  {ffo+1»7rffo+1>0?+1}  =  {°»°>°}  ’  or  ,s  small  -  take  p,l+1  •  n|  as  boundary 
control;  if  not  compute 


(3.45) 
and  then 


/n(lgo+1l2  +  No4-1!  + 


(3.46)  {u/*+1,3TwJ+1,U»J  +  1}  =  {ffo  +  1»7rgo+1.gj1+1}  +7nK,™o,<}- 
Do  n  =  n  +  1  and  go  to  (3.36). 

Remark  S.l  :  Problems  (3.33)  and  (3.38)  are  particular  cases  of 

(3.47) i  f  (u  -  *7*p)udx=  f  fvdx,\fveL2[Cl), 

J  n  Jn 

(3.47) 2  j  (p  •  q  +  uV  •  q)dx  =  0,  Vqe#(n,  div), 

J  n 

which  is  the  mixed  formulation  of  the  following  Dirichlet  problem 


(3.48) 


— Au  -f  u  =  f  in  fl,  u  =  0  on  I\ 


Observing  that  V  •  qeL2(fi),  Vqei7(f2,  div),  we  can  eliminate  u  from  (3.47)i,  (3.47)2 
to  obtain  that  p  satisfies  (if  /eL2(fl))  : 

qI  j  peH(Q,  div), 

3  4  1  Jq(V  •  PV  '  q  +  P  •  <l)dx  =  ~  In /V  •  qdz,VqeH{n,div). 

Solving  (3.49)  (in  fact  its  discrete  variants)  is  fairly  easy  and  can  be  done  by  con¬ 
jugate  gradient  algorithms  (see,  e.g.,  (9j  for  details).  Once  p  is  known,  one  obtains 
easily  u  from  (3.47);  .  Combining  the  above  algorithm  with  the  mixed  finite  element 
approximations  and  time  discretization  schemes  of  the  wave  equation  discussed  in 
Section  2  is  (almost)  straightforward;  this  issue  will  be  discussed  in  a  forthcoming 
paper. 
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3.4.3.  Numerical  experiments. 


The  mixed  finite  element  approximation  and  time  discretization  schemes  of  the 
wave  equation,  described  in  Section  2,  have  been  combined  to  algorithm  (3.30)  - 
(3.46),  to  solve  problem  (3.11)  when  Cl  =  (0,1)  x  (0,1)  and  T  —  2\/2  .  Using  the 
Fourier  series  techniques  described  in  [13]  we  have  computed  those  initial  data  u0 
and  u i  for  which  the  solution  e(=  {e„,ei})  of  (3.11)  is  given  by 

(3.50)  ea(x  1,2:2)  =  sin  7TZi  sin7ri2,  ei  =  ir\/2 e0. 

We  have  used  the  mixed  finite  element  approximations  of  Section  2,  with  k  =  1 
and  Rh  the  regular  partition  of  Cl  associated  to  the  vertices  {ih,jh}  with  0  <  itj  < 
N,  N  being  an  integer  such  that  Nh  =  1  ;  we  have  taken  N  =  16,32,64  .  The 
time  discretization  of  the  various  wave  equations  involved  in  the  calculations  was 
obtained  using  the  (conditionally  stable)  explicit  scheme  described  in  Section  2. 
Obtaining  the  (approximate)  values  of  the  control  =  p  •  n  on  ,  was  quite 
easy  since  the  values  of  the  fluxes  (i.e.  of  the  normal  components  of  ),  at  the 
element  interfaces  and  at  the  boundary  T  ,  are  the  natural  degrees  of  freedom  for  the 
functions  belonging  to  the  finite  dimensional  space  Qh  approximating  H(Cl,div). 
For  h  =  l/l6(resp.l/32, 1/64)  the  finite  dimensional  variant  of  algorithm  (3.30) 

-  (3.46)  converges  in  48  (resp.  72,  119)  iterations  (the  number  of  iterations  varies 

-  approximately  -  like  \'N  ).  These  numbers  are  much  higher  than  those  obtained 
in  r13\  where  the  space  approximation  was  achieved  by  a  conforming  finite  element 
method,  coupled  to  a  biharmonic  Tychonoff  regularization  to  eliminate  spurious 
oscillations.  On  the  other  hand,  using,  as  in  the  present  paper,  mixed 'finite  element 
approximations,  it  is  not  necessary  to  use  regularization  to  obtain  very  good  nu¬ 
merical  results,  as  shown  in  Figures  3.1  (a),  (b),  (c)  (N=6),  3.2  (a),  (b),'  (c)  (N=32), 
3-3(a;,  /b),  (c)  (N=64). 

Figures  (a)  (resp.  (b))  show  the  variation  of  the  exact  (-)  and  computed  (•  •  -)e0 
frcsp.  €1  ),  for  0  <  x\  <  1,^2  =  -5.  Figures  (c)  show  the  variation  on  (0,!T)  of  the 
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£2(r)~  norm  of  the  exact  and  approximate  boundary  controls. 

All  the  above  calculations  have  been  done  on  a  CRAY  X-MP  supercomputer. 

4.  CONCLUSION. 

In  this  paper  we  have  discussed  the  application  of  mixed  finite  element  methods 
to  the  numerical  solution  of  direct  or  inverse  problems  for  time  dependent  equations. 
These  mixed  methods  are  robust  and  accurate.  They  are  however  more  complicated 
to  implement  than  the  traditional  finite  element  methods.  Indeed  many  important 
issues  remain  concerning  the  practical  use  of  the  mixed  methods  considered  here, 
such  as  speeding  up  calculations  by  multigrid  and/or  domain  decomposition  meth¬ 
ods  (cf.  [10]);  we  intend  to  investigate  them  in  the  near  future. 
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